!************************************************************
! this subrouitne calculate invese Symmetry point of
! isentropic supersonic flow
!************************************************************
! Created by : B. G.
! Date       : 2015-08-13
! Revised    :
!************************************************************
! Nomenclature: 
!   Vp~theta0	properties along l-c-line(Mach line, 24, *p) and streamline(34, *0)
!   lp~l0      coefficients of functions, *p--l-c-line, *0-streamline
!   *4f        properties of previous loop of point 4
!   atemp	   temporary speed of sound, m/s
!   i          loop indicator
!   
!************************************************************
! Warning :
!   1. angle in rad, not degree
!************************************************************

subroutine FreeP()
   !       4()
   ! 3
   !     2
   use VariableDef
   implicit none
   ! given variables
   ! real*8::x2,y2,V2,theta2,p2,rho2,x3,y3,V3,theta3,p3,rho3
   ! temporary variables
   real*8::pp,Vp,thetap,rhop,yp,alphap,Mp,theta0
   real*8::lp,Qp,Sp,Tp,l0,atemp
   real*8::x4f,y4f,theta4f
   integer::i
   ! calculate variables
   ! real*8::x4,y4,theta4,T4,M4,Vv4,Vu4
   
   p4=p3
   V4=V3
   rho4=rho3
   
   pp=p2
   Vp=V2
   rhop=rho2
   thetap=theta2
   yp=y2
   
   theta0=theta3
   
   x4f=0.0
   y4f=0.0
   theta4f=0.0
   
   do i=1,icor
      ! determine coefficients along l-c-line, 24
      l0=tan(theta0)
      call Thermo(pp,rhop,Vp,atemp,T4,Mp)
      alphap=asin(1/Mp)
      lp=tan(thetap+alphap)
      Qp=sqrt(Mp**2-1)/(rhop*Vp**2)
      Sp=delta*sin(thetap)/(yp*Mp*cos(thetap+alphap))
      ! determine x4, y4, Tp and theta4
      x4=(-lp*x2+l0*x3+y2-y3)/(l0-lp)
      y4=y3+l0*(x4-x3)
      Tp=-Sp*(x4-x2)+Qp*p2+theta2
      theta4=Tp-Qp*p4
      ! justify if point 4 arrived
      if ((abs(x4-x4f) <= el) .and. (abs(y4-y4f) <= el) .and.&
         &(abs(theta4-theta4f) <= etheta*theta4f)) exit
      
      pp=(p2+p4)/2
      rhop=(rho2+rho4)/2
      Vp=(V2+V4)/2
      thetap=(theta2+theta4)/2
      yp=(y2+y4)/2
      theta0=(theta3+theta4)/2
      
      x4f=x4
      y4f=y4
      theta4f=theta4
      
   end do
   call Thermo(p4,rho4,V4,atemp,T4,M4)
   Vv4=V4*sin(theta4)
   Vu4=V4*cos(theta4)
   
end subroutine FreeP
